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Models for the mass function and assembly histories of dark halos: 
an approach to inventory isolated overdense regions in random fields 
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ABSTRACT 

In order to attain a statistical description of the evolution of cosmic density fluctuations in 
agreement with results from the numerical simulations, we introduce a probability conditional 
formalism (CF) based on a complete inventory of isolated overdense regions in a density 
random field. This formalism is a useful tool for describing at the same time the mass function 
(MF) of virialized dark haloes, their mass aggregation histories (MAHs) and merging rates 
(MRs). The CF focuses on virialized regions in a self-consistent way rather than in mass 
elements, and it offers an economical description for a variety of random fields. Within the 
framework of the CF, we confirm that, for a Gaussian field, it is not possible to reproduce 
at the same time the MF, MAH, and MR of haloes, both for a constant and moving barrier. 
Then, we develop an inductive method for constraining the cumulative conditional probability 
from a given halo MF description, and thus, using the CF, we calculate the halo MAHs and 
MRs. By applying this method to the MF measured in numerical simulations by Tinker et 
al. (2008), we find that a reasonable solution, justified by a mass conservation argument, is 
obtained if a rescaling -increment by ~ 30%- of the virial mass defined in simulations is 
introduced, and a (slight) deviation from Gaussianity is taken into account. Thus, both the 
MAH and MR obtained by a Monte Carlo merger tree agree now with the predictions of 
numerical simulations. We discuss on the necessity of rescaling the virial mass in simulations 
when comparing with analytical approaches on the ground of the matter not accounted as 
part of the halos and the halo mass limit due to numerical resolutions in the simulations. 
Our analysis supports the presence of a diffuse dark matter component that is not taken into 
account in the measured halo MFs inasmuch as it is not part of the collapsed structures. 

Key words: cosmology: dark matter — large-scale structure of Universe — galaxies: haloes 
— galaxies: formation — methods: statistical — methods: analytical 



1 INTRODUCTION 

According to the contemporary cosmological paradigm, cosmic 
structures emerge from the gravitational growth of primordial dark 
matter density fluctuations. A central problem in the last four 
decades has been the connection between this primordial fluctu- 
ation field and the abundance and assembly history of the virial- 
ized dark matter haloes; in more detail, their mass function (MF), 
mass aggregation histories (MAHs), and merger rates (MRs). Since 
the seminal work by Press & Schechter (1974; hereafter PS), many 
statistically-based analytical formalisms were developed in order to 
perform such a connection (see for reviews e.g., Zentner 2007; Mo, 
van den Bosch & White 2010). 

In order to achieve an association between the density fluctua- 
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tion field and the virialized haloes, one needs: a statistical descrip- 
tion of the density fluctuation field, an inventory of the overdense 
regions which will be associated to the virialized objects, and a 
model for the dynamical evolution of the overdensities, including 
an operational criterion of collapse (virialization). 

PS assumed a Gaussian density field and used the spherical 
collapse model (Gunn & Gott 1972). A region of the density field is 
assumed to end as a virialized halo when its linearly evolving over- 
density exceeds a critical value 5 C . The virialized mass at a fixed 
scale is assumed to be given by the contributions of the regions of 
this scale overdense by S c (the PS Ansatz) plus the underdense re- 
gions of the same scale contained in larger regions overdense by S c ; 
the latter statement enunciates the so-called cloud-in-cloud prob- 
lem. In order to account for all mass, PS assumed that the second 
contribution equals to the first one, justifying this the fudge factor 
of 2 introduced in their inferred halo MF. 

Based on a more rigorous statistical description, the excursion 
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set (ES), also know as extended PS (EPS) formalism (Peacock & 
Heavens 1990; Bond et al. 1991), allows to overcome the cloud- 
in-cloud problem. The ES formalism applied to a Gaussian field 
provides a tool to compute, besides the unconditional MF, the con- 
ditional MF that can be used for generating halo merger trees (Bond 
et al. 1991; Lacey & Cole 1993). In this formalism, the result de- 
pends critically on the filter used to average the linear density fluc- 
tuation field at different mass scales. The mathematical solution is 
straightforward when a sharp fc-space filter is adopted. For this fil- 
ter, the PS result including the fudge factor of 2 is recovered. In 
the calculation of the merger trees, the sharp fc-space filter implies 
independent steps (Markovian random walks) along the mass tra- 
jectory. Note that the ES formalism has a conceptual problem for 
predicting the halo mass in which a particular mass element ends 
up (Mo et al. 2010, §§7.2.2), problem that may lead to an inaccurate 
buildup of MAHs and MRs. 

With the advent of large cosmological N-body simulations, 
the whole non-linear process of dark matter gravitational evolution 
and collapse into virialized haloes could be followed, though with 
strong limitations due to mass resolution. Do approaches based 
on the Gaussian ES formalism allow to describe correctly the re- 
sults from the simulations? In particular, do results from these ap- 
proaches agree at the same time with the MF, MAHs, and MRs of 
haloes as measured in simulations? 

A non-negligible discrepancy between the MFs obtained in 
the ES formalism and the N-body simulations has been early re- 
ported. The introduction in the ES formalism of the elliptical grav- 
itational collapse instead of the spherical one (mass-dependent in- 
stead of constant 5 C , respectively) helped to overcome this problem 
(Sheth, Mo & Tormen 2001). However, in this case the ES formal- 
ism does not provide an analytic formula for the conditional MF, 
and the merger trees based on the Gaussian ES formalism show 
deviations with respect to simulations in the progenitor mass distri- 
butions (Sheth & Tormen 2002), the mass contained within all the 
progenitors (Neistein et al. 2006), and the average main progeni- 
tor MAHs and MRs (Wechsler et al. 2002; van den Bosch 2002). 
On the other hand, direct measures of the conditional MFs in N- 
body simulations show that they depart from the corresponding 
functions calculated with the ES formalism for a Gaussian field 
(Cole et al. 2008; Neistein et al. 2010). Thus, a possible source of 
the discrepancies lies in the strict assumption of Gaussianity in the 
ES formalism. Note that the question is not about primordial non- 
Gaussianity; the N-body simulations use indeed a Gaussian density 
field that is evolved analytically (e.g., by means of the Zel'dovich 
approximation) until the quasi-linear regime is reached. The point 
is that non-negligible deviations from Gaussianity at the scales of 
interest have likely happen already during this regime (Coles & 
Jones 1991). Here, we set aside Gaussianity and the ES formalism 
and handle the problem of the inventory of overdense regions in 
arbitrary random fields. 

We present an approach aimed to state the problem of the 
inventory of overdense regions in random (Gaussian or non- 
Gaussian) fields from a very general point of view. We assume 
that the density fluctuation field is fully characterized by its field 
conditional probability function. In order to face the problem of as- 
sociating the halo mass to the overdensity in an alternative way to 
the ES formalism (see the reference Mo et al. 2010 cited above) we 
proceed as follows. In our conditional formalism (hereafer CF), we 
make use of the concept of isolated regions introduced by Jedamzik 
(1995; see also Yano, Nagashima & Gouda 1996; Nagashima 2001) 
as a natural way to connect overdensities to virialized regions. 
We enunciate the isolated overdense regions inventory theorem, 



through which the given field conditional probability function of 
the random field is linked to the inventory of isolated overdense re- 
gions inside larger isolated regions of lower overdensity, being the 
former those that eventually collapse into virialized haloes if their 
overdensities equal to 8 C . 

The CF offers an alternative EPS formulation. It gives an easy 
way to develop a Monte Carlo algorithm for calculating the MAH 
and MR of the growing virialized haloes, because it supplies natu- 
rally the progenitor conditional probability function \j of finding a 
specific isolated region (progenitor) inside a given isolated region 
(descendent). When the problem of finding isolated overdense re- 
gions inside larger ones is extended to the overall Universe, then 
such conditional probability reduces to an unconditional probabil- 
ity, which gives the halo MF. Nevertheless, if one starts from an 
unconditional probability function, for example on the basis of a 
halo MF obtained from N-body simulations, and tries to go back to 
the conditional probability function, in order to build merger trees 
compatible with the given MF, then multiple possibilities appear. 
The CF offers useful tools to make a choice. As we will discuss in 
§§4.2, we decide to renounce to the Gaussian hypothesis even if, 
for simplicity, we retain the (Gaussian) PS Ansatz in a generalized 
version. 

The CF might be a powerful tool for economically describing 
in a consistent way the simulation results related to the MF, MAHs, 
and MRs of the virialized haloes. These descriptions can then be 
easily implemented in semi-analytical and semi-empirical models 
of galaxy evolution, and can be extended to masses and epochs, 
where the resolution is a limit for the simulations. On the other 
hand, the CF may allow to explore economically cosmic structure 
formation in alternative cosmologies or in cases where the density 
fluctuation field is intrinsically non-Gaussian. Our approach could 
be applied also in other astrophysical problems, where random den- 
sity fields are introduced, e.g., the formation and evolution of dense 
gas structures in the interstellar medium as giant molecular clouds, 
massive clumps and cores (Hopkins 2012). 

This paper is organized as follows. Section 2 is devoted to 
present the CF; in §§2.1, an equation for the isolated overdense 
regions inventory is derived, and in §§2.2 the algorithms to build 
up the merger trees are introduced. Section 3 deals with the par- 
ticular case of Gaussian random fields. The halo MF is derived in 
§§3.1, while the MAH and MR are calculated in §§3.2. A discus- 
sion on the moving barrier effects is presented in §§3.3. In view 
of the discrepancies obtained when using as input the Gaussian 
statistics, in Section 4 an heuristic approach, based on the use of 
the halo MF from simulations as input, is presented. In §§4.1 we 
approach the problem of the difference on the mass estimates in 
the analytic method and in the simulations, and remark the ne- 
cessity to introduce a mass rescaling, representative of a diffuse 
matter component, to connect the two standpoints. The results con- 
cerning MAHs and MRs are presented in §§4.2. Our conclusions 
are given in Section 5. We adopt a flat ACDM cosmology with 
fi M = 0.27, Qa = 0.73, h = 0.71, a s = 0.9, n = 1. 



2 THE CONDITIONAL FORMALISM 

The ES formalism provides a means to derive the halo MF, account- 
ing for the cloud-in-cloud problem, and to build up the halo MAHs. 



1 For brevity, hereafter we will omit the specification "progenitor" when 
referring to the conditional probability. 
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However, as mentioned above, the Gaussian ES formalism seems to 
have difficulties for predicting at the same time the MF, MAH, and 
MRs in agreement with numerical simulations. Also, with the ES 
formalisms it is not easy to deal with non-Gaussian random fields 
(Inoue & Nagashima 2002), though some progress has been made 
recently by introducing a stochastic barrier and non-Markovian cor- 
rections (Maggiore & Riotto 2009), by manipulating the step-size 
distribution of the random walks (Lam & Sheth 2009) or by ac- 
counting for correlations between steps (Musso & Paranjape 2012). 
On the other hand, being the ES formalism focussed on a generic 
mass element, the virialized halo mass is associated to a smooth- 
ing scale, which may lack of a real physical meaning. In order to 
overcome these shortcomings, we introduce below the conditional 
formalism (CF) inspired by the Jedamzik (1995) approach. This 
formalism is based on a complete inventory of the isolated over- 
dense regions of a random density field described by a cumulative 
field conditional probability function. 



Following, a theorem focussed on an isolated overdense re- 
gions inventory and aimed to establish a connection between the 
cumulative field conditional probability, F(M' , S'\M, 8), and the 
conditional MF, N(M' , S'\M, 8) with 8' > 8, is presented. 

Consider a region S with overdensity 8 and mass M, De- 
fine an arbitrary mass M[ ^ M. A decomposition of the mass 
range (M[, M) into n ordered steps, dM[, identifies a sequence of 
dNi = J\f(Mi,8'\M, 5)dM'i isolated regions 1(5') with overden- 
sity 8' and masses between M[ and M[ + dM- inside S. 

For any given 8" ^ S', the amount of mass in S assembled by 
regions with mass M[ and overdensity between 8" and 8" + dS" is 

M<j>(M{,8"\M,d)dS" 

This mass is also given by the sum of the regions of mass M[ with 
overdensity between 8" and 8" + dS" located in each one of the 
isolated regions 1(8') with masses between M[ and M[ + dM[, 
and overdensity 8' of the above defined sequence: 



2.1 A complete inventory of isolated overdense regions 

Let start with some key definitions. 

Isolated regions: an isolated region with overdensity 8 is a con- 
nected region not included in a larger connected region with over- 
density ^ 8, while at the same time, outside of the totality of iso- 
lated regions with overdensity 8 there does not exist any region with 
overdensity ^ 8. Hereafter we will refer to an isolated region with 
overdensity 8 by means of the symbol I(<5)q 
Virialized regions: a virialized region is an isolated region whose 
overdensity 8 is equal to a critical value S c . 

4>(M', S'\M, S)dS': is the field conditional probability to find a re- 
gion of mass M' with an overdensity between 8' and 8' + dS', 
contained inside a larger region of mass M and overdensity 8. 
J\[(M',8'\M,8)dM'; is the number of isolated regions, i.e. the 
conditional MF, with overdensity 8' and mass between M' and 
M' + dM', contained inside a larger region S with overdensity 8 
and mass M. 

The cumulative field conditional probability to find a region of 
mass M' with overdensity ^ 8' inside a larger region of mass M 
and overdensity 8 is 

F(M',S'\M,S) = / 4>(M',S"\M,S)dS" 

JS' 

Our CF concerns random fields that can be defined by such condi- 
tional probability. 



2 The definition of isolated region allows to introduce new mathematical 
concepts. Given a space with a distribution of isolated regions I n (<5), one 
can cut an arbitrary region of such space along a specific border. The con- 
sequence will be that some I n (<5) will be internal, others will be external 
and finally some of them will be on the border of that region. By defini- 
tion, when no one of the isolated regions l n (8) will be on the border, such 
region will be named suitable with respect to its internal isolated regions 
overdense by 8. It is straightforward to see that an isolated region 1(5) is 
a suitable region with respect to all its internal isolated regions more over- 
dense than <5. Another example is given by an isolated region 1(5) where 
one (or more) internal isolated region(s) I r (<5') with 8' > S has (have) 
been removed. The result is a suitable region S with respect to its internal 
isolated regions 1(8'), its overdensity is less than 8, and typically it will not 
be an isolated region. Any region of S more overdense than <5' is contained 
inside an internal isolated region 1(8'). Afterwards when we will refer to 
a region S actually we will refer to a region suitable with respect to the 
internal isolated regions overdense by 8 involved by the problem. 



J2 M-<j)(M[,S"\M-, S')dS" dNi 

No other contribution exists because outside the totality of the iso- 
lated regions 1(8') in S there are not regions with overdensity ^ 8' . 
By equating the last two contributions, we can write 

A/f 
0(M{,8"\M,8)d8" = Y^-j^ d -Vi <f>(MiS"\M-,8')d8", 

i=i 
from which we derive the following integral equation: 

<j>(M',S"\M,8)dS" = 



w, Af " 



dM" -j- N(M" ', 8'\M, 8) 4>(M', S"\M", 8')d8" 



By integration on 8" between 8' and oo, the previous equation 
writes in terms of the cumulative field conditional probability: 

F(M',8'\M,8)= (1) 

dM ' ~W F( - M '' 5 '' M "' <5 '- ) A/ '( M "' 5 '\ M ' 5 ) 

An alternative formulation of Eq. ([T} is obtained introducing 
the conditional probability p(M' , 8'\M, 8)dM' of finding an iso- 
lated region of mass between M ' and M' + dM' and overdensity 
8' inside a region S. Given a region S with overdensity 8 and mass 
M, our interest is now to identify inside it an isolated region with 
overdensity 8' . The quantity 

AM" = M"M(M", 8'\M, 8) dM" 

gives the amount of mass of S assembled by isolated regions with 
overdensity 8' and mass between M" and M" + dM" . There- 
fore, the probability that an arbitrary mass element of S is part of 
an isolated region with overdensity 8' and mass between M" and 

M" + dM" is 



p(M",8'\M,d)dM" = 



AM" 
M ' 



that is 



p(M",8'\M,S) 



AT 
M 



N(M",8'\M,S). 



(2) 



The previous statement is equivalent to say that 
p(M",S'\M, 8) dM" is the probability to find an isolated 
region with overdensity 8' and mass between M" and M" + dM" 
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inside S. Introducing Eq. ((2} into Eq. (QJ, we obtain 
F(M',5'\M,S) = 



dM" F{M', S'\M", 8') p(M", S'\M, 8) 



(3) 



which is basically a Volterra integral equation. 

The fundamental aspect of Eq. Q} is that for any given cu- 
mulative field conditional probability F(M' , 8'\M, 5), it gives the 
conditional MF, Af(M',6'\M, 5). Its main difference with the 
Jedamzik (1995) result is that Eq. l[T} is applied to any region S 
with overdensity S and mass M and not only to the overall uni- 
verse. For this reason the cumulative field conditional probability 
F(M' , S'\M, S) appears now on both sides of the equation. Such 
property highlights the intrinsic meaning of the CF It is interesting 
to remark that the previous isolated overdense regions inventory 
theorem solves automatically the cloud-in-cloud problem. The Eq. 
(0 represents the generalization of the PS formalism. Such formal- 
ism may be recovered for the special random fields wherein 



F(M',S\M,S) = k 
with k constant. In such case Eq. Q} writes 



p{M',S'\M,S) 



1 dF(M',5'\M,8) 
"k dM' 



(4) 



(5) 



The PS Ansatz in terms of the conditional probability and including 
already the fudge factor of 2, is obtained in the particular case of 
k = 1/2 (corresponding namely to Gaussian fields, see below). 

A random field for which Eq. (|4j is satisfied is very simple to 
handle in terms of Eq. l|5}. Note that not only Gaussian fields obey 
Eq. $4$ but also a variety of non-Gaussian fields. Therefore, in order 
to find an analytical formalism in agreement with simulations and 
to avoid the complex numerical integration of Eq. (Q3, it is natural 
for us to look at first among random fields which satisfy Eq. l(4). 
This is the strategy followed in Section 4. 

2.2 Merger trees build up 

A merger tree realization needs to identify progenitors of a given 
descendant. This is done by means of a Monte Carlo algorithm that 
uses the conditional probability and MF. The choice of a Monte 
Carlo algorithm is not unique and may hide subtle physical and 
technical questions (for a review see Zhang et. al. 2008). The sim- 
plest criterion establishes that the Monte Carlo algorithm has to 
reproduce the conditional probability and MF predicted by the the- 
ory. However, this criterion is not sufficient, because the overall 
meaning of the involved probability, i.e. the stocahstic nature of the 
problem, has to be taken in to account. We will introduce a Monte 
Carlo algorithm optimized to build up a merger tree. Afterwards 
we will check our algorithm with tests that involve the conditional 
probability and MF, as well as the MAH and MR in order to verify 
also the agreement with the stochastic nature of the problem. 

Only in this section, to ease the reading of the formulae, we 
will indicate the variables of the extracted regions by a lower case 
letter. 

Integrating p(m' ,8'\M, S)dm' between and m and identi- 
fying the obtained cumulative probability with a random number a 
uniformly distributed in the interval [0-1], we write 



p(m' , S'\M, 8) dm' , 



(6) 



which represents the basis for a Monte Carlo approach to build up 
merger trees. 



Given an isolated region Io(5o) with overdensity So and mass 
Mo, the extraction of the random number a\ by Eq. ® allows to 
identify inside it a first isolated region ii(5') of mass mi with a 
given overdensity 8' . Taking into account that the isolated over- 
dense regions inventory theorem has been proved for general re- 
gions (see footnote 2 ), after the identification of the isolated region 
ii(<5'), the extraction may continue on the complementary region 
Si(5i) with mass M\ — Mo — m\ and volume Vi = Vo — Vi 
obtained from Io(5o) removing ii(<5') from its inside. In terms of 
the density, from the volume relation we write 



Mo 
Pa 



m i 

V 



Mi 
Pi : 



Being the region density p = p(8 + 1), where p is the background 
density and 8 the overdensity, we have 



Mo 



and finally 



So 



5i = 



7711 Ml 



Mo — mi 



Mo mi 

So + 1 ~ S' + l 

At this point the extraction of the next isolated region can be 
earned out on the region Si (Si ) with overdensity Si and mass Mi . 
To perform such extraction, for simplicity, we will use Eq. (|6]l with 
a new random number «2, applied on the region Si(<5i). By pro- 
ceeding in this way we implicitly assume that the same statistics 
of Io(<5o) holds for Si(<5i), excluding any sort of correlation. For 
a discussion of possible correlations, see Sheth & Lemson (1999). 
According to what has been said above, afterwards we will proceed 
to test the robustness of such hypothesis. 

A recurrent algorithm may be particularly useful now. From 
the n th complementary region S n of mass M n with overdensity 
S n calculated by the previous method, using Eq. ((6} with a random 
number a n +i we can identify inside it an isolated region i n +i(5') 
with mass m n +i and overdensity 8' . At this point the next comple- 
mentary region S n +i with mass 



M n+ i = M n - m n+ i 



will have the overdensity 

8n+l = 



M n - m n+1 



M n 



m n +i 



(7) 



(8) 



S n + 1 S' + l 

Equations 10 and {8]( guarantee a mass conservation and make the 
procedure particularly handy because they do not make explicit ref- 
erence to the previous i-steps. 

Equations l[6}{8} are particularly useful for constructing the 
halo merger trees. Suppose a simple case when the isolated viri- 
alized regions are defined by a critical overdensity S c function of z 
only. Given a descendant at z with overdensity S c (z), the progeni- 
tors at z + Az have overdensity 8 c (z + Az). The finding of such 
progenitors (branches) is performed by applying the Monte Carlo 
extractions of Eq. (f6]l to the (complementary) region described by 
Eqs. (|7) and ([8}. This procedure identifies a sequence of progen- 
itors at 2 + Az that can be sorted in order of decreasing mass. 
The largest (primary) progenitor identifies the main progenitor by 
which the mass aggregation history, MAH, is obtained. The overall 
population of less massive progenitors fulfills the merger tree and 
allows to find the merger rate, MR, along the MAH. 

Note that the condition p = implies 5 = — 1. Therefore, for 
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the complementary regions the progenitor extractions can be per- 
formed meanwhile S n ^ — 1. This means that the progenitors are 
extracted until all matter in the volume is exhausted, which implies 
that the entire mass contained in the descendant should come from 
virialized structures. However, the physics of halo mass assembly 
is more complex. A fraction of the mass can be present as diffuse 
matter. The strict limit 5^—1 implies lack of diffuse matter. In 
this paper we explore solutions with 8^—1, however we will have 
in mind the physical consequences of such a hypothesis. 

The above described method has been extensively used in 
the case of Gaussian CDM density fields for generating the ini- 
tial conditions (the MAHs) for the virialization process of dark 
haloes (Avila-Reese, Firmani & Hernandez 1998), and the subse- 
quent formation and evolution of disc galaxies inside them (Firmani 
& Avila-Reese 2000,2009; Avila-Reese & Firmani 2000). The ob- 
tained halo properties agree with those found in cosmological nu- 
merical simulations (Avila-Reese et al. 1999). 

It is interesting to remark that the mass distribution of the first 
progenitor obtained with the Monte Carlo method can be compared 
with the conditional probability function, while the mass distribu- 
tion of all the progenitors can be compared with the conditional 
MF This offers a first test of our Monte Carlo algorithm by com- 
paring its one-step results with the straight analytic formulae from 
the statistics (in §§3.2 we will show this point in some detail). 

It is important to highlight that the mean halo MAH has to be a 
convergence solution of the Monte Carlo algorithm for the integra- 
tion redshift step Az — > 0, otherwise it is physically meaningless. 
Because of the stochastic nature of the problem, this property is de- 
pendent on the behavior of the given cumulative field conditional 
probability. We will focus on well behaved cumulative field condi- 
tional probabilities, such that guaranty Az — > convergence in 
the MAH buildup. 



3 GAUSSIAN RANDOM FIELDS 

In the case of a Gaussian random field, the cumulative field con- 
ditional probability to find a region with mass M' and overdensity 
^ 8' contained inside the region of mass M with overdensity 8 is 
given (see Bower 1991) by: 



F(M',S'\M,S) = -j= I 



where 



7 = 



(5' - S) 



\] 2 i?w ~ ct m) ' 



(9) 



(10) 



and <xm is the mass variance calculated from the density fluctuation 
power spectrum of the random field. Making the limit 8' — > 5, Eq. 
(|9j gives 

F(M', S\M, 5) = - (M' < M), (11) 

which is Eq. © with k = 1/2. Then Eq. $5$ reduces to 

p(M\5'\M,5) = -2 dF ( M '> 5 '} M > S \ (12) 

dM' 

This is just the PS Ansatz in terms of the conditional probability 
with the fudge factor of 2 included. Introducing Eq. l(9) into Eq. 
l !12t we obtain 



P (M',S'|M,*) = ^e-^. 




Log(M/M ) 

Figure 1. Halo MFs plotted as (M 2 /p)dN/dM from Tinker08 (solid 
black lines) compared to the PS MFs obtained with the CF (solid red lines) 
for z = 0, 1, 2, 3, 4. The dotted red lines show the PS MFs with rescaled 
critical overdensity and lowered amplitude (see text). 



It is easy to verify that such conditional probability is normalized, 
i.e., its integral in M' from to M (7 from to 00) is 1. 



3.1 Unconditional mass function for haloes 

Eq. d 1 3 fc is applied to a region of mass M and overdensity 8. We can 
extend such region to the entire universe making M —> 00, o~m — > 
0, and 8 — > 0. Hereafter, in such limit case of the entire universe, 
the index ' will be omitted. A recipe for the collapse criterion may 
be obtained from the top-hat spherical gravitational collapse, which 
establishes that an isolated region at redshift z is virialized when its 
overdensity is equal to a critical value 5 C (Navarro Frenk & White 
1997). For a flat cosmology with cosmological constant (Qm + 
Oa = 1): 



S c = 



1.686 n° M 0055 



(14) 



D+(z) 
where D + (z) is the linear growth factor normalized to 1 at z = 



By defining the number density per mass unit 

dN p , , , - , „> 
— = i^(MA|oo,0), 

definition that is not limited to Gaussian fields, and introducing 

_tfo_ 

CM ' 

we obtain the halo MF 



(15) 



(16) 



dv 



2 _Z. 

e 2 

n dM 



(17) 



(13) 



M dN_ 
p dM 

which coincides with the PS MF. 

In Fig. [T] the solid red curves show the Gaussian-case (PS) 
halo MF obtained from Eq. {FT) at 2 = 0, 1, 2, 3, 4, while black 
solid curves are the accurate fitting formulae to cosmological N- 
body simulations provided by Tinker et al. (2008; hereafter Tin- 
ker08); the parameters for the virial mass at the overdensity A(z) 
(Bryan & Norman 1998) corresponding to our cosmology were 
used. The comparison reveals the well known excess (deficit) of 
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intermediate (high) mass haloes predicted by the PS formalism. 
The Tinker08 fitting curves are valid for halo masses above ~ 
10 h~ Mq, which corresponds to the minimal mass at which 
haloes are reasonable resolved in the simulations studied by these 
authors; this limitation makes uncertain the normalization of the 
overall MR Just in order to explore possibilities to approximate 
the PS MF to the simulations MF, 8 C has been rescaled by a factor 
0.86 to fit the high mass cut-off of the N-body (Tinker08) curves 
(for similar earlier attempts see e.g., Carlberg & Couchman 1989; 
Klypin & Rhee 1994), and the normalization has been lowered 
by a factor 0.61 to fit the maxima (red dotted curves). Note that 
the PS MF is normalized independently of S c , and that the dot- 
ted curves imply now a not normalized MF. It is evident that in 
spite of these mass-independent transformations of the PS MF, the 
comparison with simulation results continue failing. The exercise 
presented above is congruent with the one carried out in Sheth & 
Tormen (1999), where they showed that besides these operations, 
the PS MF should be multiplied by a v (mass)-dependent factor in 
order to fit the numerical simulations analyzed by them. 

3.2 Mass aggregation histories and merger rates 

For a Gaussian field and the critical overdensity given by Eq. \\A\ , 
if M' -> 0, then a M , ->• oo, 7 ->• and F{M',6'\M,6) -> 1/2. 
Using Eqs. il_t and ([9), Eq. ® writes as 



— / e~ ? d£. 



(18) 



This operation is equivalent to make a random choice of a number 
7 with a normal deviate, zero mean and variance 0.5. Hence, using 
Eq. d lOt and making reference to §§2.2, we obtain: 



2 



(Si - Si) 2 

2j 2 



2 



(19) 



where the symbols are the same as in §§2.2, and So = 5 c (z) is the 
(rescaled critical) overdensity of the descendant at redshift z, while 
8[ = S c (z + Az) is the (rescaled critical) overdensity of a progen- 
itor at redshift z + Az. The i sequence defines the complementary 
regions of mass Mi with overdensity Si, as well as the progeni- 
tors of mass Mf with overdensity 8[. Equations Q, $fQ and <19t , 
togheter, define the Monte Carlo algorithm to build the merger tree 
of a virialized halo with a given mass at a given redshift. 

In Fig. [2] we plot the average of 2 10 4 different MAHs cor- 
responding to a present-day halo of mass Mh = 10 13 Mq (red 
solid line; in this case, the critical overdensity is not rescaled). As 
in Fig.Q] the red dotted line corresponds to the rescaling in the crit- 
ical overdenstiy and the amplitude-reduction correction of the MF 
mentioned above. The dashed line represents the fit to the corre- 
sponding average MAH measured in the Millennium Simulations 
(Fakhouri, Ma & Boylan-Kolchin 2010). 

For present-day haloes of Mh = 10 13 Mq, Fig.[3]shows the 
mean MR histories per unit of redshift with merger ratios greater 
than £ = M s /M p (M p and M s are the masses of the primary 
and secondary progenitors, respectively). From bottom to top, £ > 
0.3 and 0.03 (green lines) and £ > 0.1, 0.01, and 0.001 (blue 
lines). The coding of solid and dotted lines is the same as in Fig. 
[2] The average halo MRs as a function of the £ threshold at z = 
compare well with the results reported in Fig. 3a by Fakhouri et 
al. (2010). The small change of these mean MRs with z is also in 
general agreement with the numerical simulation results reported 
in Fakhouri & Ma (2008) and Fakhouri et al. (2010). 

The main conclusion from the exercise presented here is that 




Log(1+z) 

Figure 2. Average MAH (solid red line) for a Mo = 1O 13 M0 present- 
day halo for the PS case, and the same for the amplitude-reduced case with 
rescaled critical overdensity (dotted red line, see text). The dashed black 
line shows the corresponding mean MAH obtained from the Millennium 
Simulations (Fakhouri et al. 2010). 
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Figure 3. Average MRs for the PS case (from bottom to up, green lines 
for £ > 0.3 and > 0.03, blue lines for £ > 0.1, 0.01 and 0.001) along 
the average MAHs of Fig. ff] (solid lines for the main case and dotted lines 
for the amplitude-reduced case with rescaled critical overdensity, see text). 
Taking into account the mass variation, this result compares well with the 
Millennium Simulations at z Ri (the thick in the vertical axis indicate the 
corresponding z = MRs reported in Fakhouri et al. 2010). The evolution- 
ary behavior also agrees with these simulations. 



for a Gaussian density field (PS), the predicted average halo MAH 
and MRs agree reasonably well with the results from numerical 
simulationsQ However, the predicted halo MF, as it is well known, 



3 For the amplitude-reduced case with rescaled critical overdensity, the 
MAH and MRs do not change significantly (Figs. [2] and [3j- In what fol- 
lows, unless otherwise stated, we will use the rescaled critical overdensity 
case. 
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Figure 4. The graphic shows the conditional probability p(M' ,8'\M,S) 
(red, left abscissa) and the conditional mass function Af(M',6'\M,S) 
(green, right abscissa) per unitary progenitor mass natural logarithm as a 
function of the progenitor mass. The dots show the average distribution ob- 
tained by the Monte Carlo algorithm for the first progenitor (bottom) and for 
the entire progenitor collection (top) using Eqs. 10 and (|8j. The descendant 
mass is 1O 13 M0, the redshift z = while the redshift step Az = 0.05 
(upper panel) and Az = 0.005 (lower panel). 
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Figure 5. The graphic shows the conditional mass function 
Af(M', S'\M, 6) (green) per unitary progenitor mass natural loga- 
rithm as a function of the progenitor mass. The dots show the average 
distribution obtained by the Monte Carlo algorithm for the entire progenitor 
collection using Eqs. Q an d GO- m me upper panel the descendant masses 
are 1O 11 M and 10 15 M Q , the redshift z = and the redshift step 
Az = 0.05, while in the lower panel the descendant mass is 10 13 Mq, the 
redshift z = 10 and the redshift step Az = 0.05. 



has an excess at intermediate masses and a deficit at high masses as 
compared to simulations at the mass range they are able to resolve 
(e.g., Tinker08). 

Particularly interesting are the one-step results shown in Fig. [4] 
for a descendant mass of 10 13 Mq at z — and redshift intervals of 
Az = 0.05 (upper panel) and Az — 0.005 (lower panel), respec- 
tively. The red curve (left abscissa) shows the analytic conditional 
probability p(M', 5'\M, 5) per unitary progenitor mass natural log- 
arithm, while the green curve (right abscissa) shows the conditional 
MF N(M' , S'\M, 5) per unitary progenitor mass natural logarithm 
as a function of the progenitor mass. Here S correspond to z — 
while S' to 2 = Az. The dots show an average distribution obtained 
with our Monte Carlo algorithm for the first progenitor (bottom re- 
gion) and for the entire progenitor collection (top region) . The total 
number of extractions, not reported here because unnecessary, is 
regulated in each case to reduce the scatter at reasonable levels. 
The marginal differences appearing in the graphics can arise from 
numerical effects, as well as from some departure from the hypoth- 
esis made in §§2.2 about the complementary regions. Fig.[5]shows 
a similar test carried out for the masses 1O 11 M0 and 1O 15 M0 at 
z — (upper panel) and for 1O 13 M0 but at z — 10 (lower panel), 
being in both cases Az = 0.05. We conclude that the agreement 
between the analytic conditional probability and MF and the Monte 
Carlo results as well as the comparison of the MAHs and MRs with 
simulations (see above) are fully satisfactory, and it indicates that 
our merger tree algorithm works properly. For our main calcula- 
tions we will assume a redshift step close to 0.05, while for step 
invariance tests we will assume 0.005. 



3.3 The effect of mass dependence in S c 

So far, we have used the spherical collapse critical overdensity 
(Eq. I14t . which is function of z but not of mass; this corresponds 
to a. fixed barrier in the ES formalism. The adoption of a mass- 
dependent critical overdensity, i.e. a moving barrier, compatible 
with an ellipsoidal collapse, improves the agreement between the 
analytic and the numerical results regarding the halo (uncondi- 
tional) MF (Sheth et al. 2001; Sheth & Tormen 2002). We extend 
now the test to the halo MAH and MR. We assume the same Gaus- 
sian random field as in Section 3 described by Eqs. d9i-(113t and 
Eqs. J15M17K and adopt the moving barrier 



S ec = a5 c (l + /(x)) 



(20) 



where \ = &m/5c, f is a function that in Sheth et al. (2001) is 
/ = bx c , and a, b, and c are free parameters. The CF can be used 
now to explore the effects of a mass-dependent critical overdensity 
on the MF, MAH and MR, under the Anstaz that the cumulative 
field conditional probability is the same as given in eq. ((9). Using 
Eq. l |17t with the parameter values given in Sheth et al. (2001), 
a = 0.87, b — 0.47 and c = 1.23, we recover almost the same 
halo MF proposed by Sheth & Tormen (1999) as a fit to numerical 
simulations, and justified in Sheth et al. (2001) as an indication 
of elliptical collapse (Fig. [6] note that both the obtained MF and 
the Sheth & Tormen one are close to the MF given in Tinker08). 
The excellent agreement shows that our Ansatz is quite reasonable. 
Within the framework of the ES approach, recent work on random 
walks with correlated steps suggests that substituting 5 ec into Eq. 
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Figure 6. Halo MFs plotted as (M 2 /p)dN/dM from Tinker08 (solid 
black lines) compared with the MFs obtained with the CF for <5 C depen- 
dent on mass (dashed red curves) for z = 0, 1, 2, 3, 4. The dotted red lines 
show the Sheth-Tormen (1999) MFs. 



dlOl l could be indeed a good approximation (c.f. Paranjape, Lam & 
Sheth 2012; Musso & Sheth 2012). 

To find the MAH and MR, we use Eqs. f9b— fT2t with 8 ec given 
by Eq. {20]». For M' ->■ (cr M > -*■ oo), 7 ->• only if c < 1, while 
7 — > ab/\/2 when c = 1 and 7 — > 00 for c > 1. In these three 
cases, ip = F(0,5' ec \M, 8 ec ) is 1/2, a value between and 1/2, 
and 0, respectively. Using Eq. l!12t . Eq. $6^ now writes 

2 



a - 2<p + 1 = 



(% 



(21) 



This equation does not have a solution when c > 1, while it does 
have a solution when c < 1 for each value of < a < 1, and 
for c — 1 only if 2ip — 1 < a < 2ip. Within the framework of 
our approach, we conclude that a mass-dependent critical overden- 
sity, as the one inferred from the moving barrier analysis suggested 
in Sheth et. al (2001), does not offer a consistent description of the 
halo MAH and MR, though, the limiting "square-root" moving bar- 
rier (c = 1) suggested by Moreno et al. (2008) in principle could 
offer a solution. 

The results obtained by the Monte Carlo method using the 
square-wot moving barrier show that both the average halo MAH 
and MR depart significantly from those found in the simulations. 
The MAH is now too extended towards the past, revealing a low 
contribution of major mergers in the halo assembly. The major MRs 
(£ > 0.1) are indeed rare with respect to the simulation results. We 
have explored the obtained effects varying the a, b, c parameters 
without any interesting result. 

Nonetheless, the most serious difficulty in our formalism 
when using a mass-dependent critical overdensity is that the con- 
vergence of the Monte Carlo algorithm, when redshift step Az — > 
0, does not accomplish anymore. It is not surprising that some 
rough agreement with the results from simulations is obtained only 
when Az > 0.5, owing to the "correct" behavior of 7 as a function 
of the progenitor mass in this case (see below). Indeed, by using 
the large Az > 0.5 time step, Moreno et al, (2008) have showed 
that the moving barrier predictions agree with the simulation con- 
ditional probability. 

Simple considerations help us to understand the failure of the 
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Figure 7. The plot shows 7 as function of the mass in the case of a fixed 
barrier (dashed lines) and a moving barrier (solid lines). The descendant 
mass is 10 13 M Q , Az = 0.005, 0.05, 0.5 from bottom to top 



mass-dependent critical overdensity in our scheme. Figure|7]allows 
to understand the reduction of the major merger rate. For a fixed 
barrier (8 C independent of mass; dashed lines), 7 increases mono- 
tonically with mass and Eq. d 1 3 fc leads to a conditional probability 
that makes possible the existence of low-mass progenitors. For the 
"square-root" moving barrier (Moreno et al. 2008; solid lines), 7 
shows a decreasing shape on a very extended low-mass range of 
progenitors; here the conditional probability Eq. (1131 is physically 
meaningless and no low mass progenitors are possible. Only in a 
small range immediately below the descendant mass, 7 increases 
with mass, allowing at this range for a physical conditional proba- 
bility, such that the progenitors are more massive than in the fixed 
barrier case. When Az > 0.5 the "correct" monotonic increasing 
of 7 with the progenitor mass is recovered. Furthermore, a neces- 
sary condition to obtain a redshift-step convergence is that 7 — > 
when Az — > with M' < M. This tendency is not fulfilled for the 
moving barrier case. 

We conclude that, when using a mass-dependent critical over- 
density (moving barrier) in our approach, the halo MF is correctly 
described but the predicted halo MAH and MR fail in reproducing 
the simulation results. An exhaustive analysis of this shortcoming 
and the reason behind it is deserved for a future work. 



4 A SEMI-EMPIRICAL STRATEGY 

In view of the difficulty of proposing a random fluctuation field 
described by an analytic cumulative field conditional probability 
function and able to generate at the same time the halo MF, MAH, 
and MR obtained in the numerical simulations, we turn on to an 
inductive (semi-empirical) approach: the probability function will 
be inferred numerically from the Tinker08 halo MF, which was 
obtained from simulations. We split the problem into two parts. 
Firstly, in §4.1 we extend the unconditional probability correspond- 
ing to the Tinker08 MF (defined only above ~ 10 11 Mq) down to 
small masses by an algebraic extrapolation, taking care that such 
an extrapolation does not imply a sharp change in the functionality; 
furthermore, we introduce a mass rescaling in order to take into ac- 
count any halo mass defect due to some diffuse matter (see §§2.2). 
Secondly, in §4.2, based on the properties of Eq. (Q], we extend the 
definition of the unconditional probability to the conditional prob- 
ability, and using it within the context of the CF, we calculate the 
halo MAH and MRs. 
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4.1 Halo mass rescaling in the numerical simulations 

Tinker08 provided an analytical fit to the halo MFs measured in 
numerical simulations at different redshifts and above a minimum 
(resolution) mass. We can calculate from these MFs the cumulative 
unconditional probability functions expressed in terms of the com- 
monly used scaled variable v or, what is the same, 7„ = v/y/2 = 
5 c {z)/V2(tm, where 7„ is the same as in Eq. d 10b but applied to 
the entire universe, i.e. for M — > oo (om — > 0) and S — > 0. 
Due to the lower limit in mass in the Tinker08 MFs, the cumulative 
unconditional probability should be constructed by integrating the 
MF from each mass M at a given redshift, corresponding to the 
argument y u , to infinity. We define the following function: 

P , x _ , r M dN 

rM(y u ;z) 

p(M,S c (z)\oo,0)dM 



— -JlT7 dM = 
dM 



(22) 



/ 

Jo 



where dN/dM, given by Eq. dl5b . is the halo MF, and the mass 
M(-y u ; z) is an intrinsic function of ~y u for a given z. As long as the 
integral in Eq. d22t is correctly normalized, P(7 U ; z) represents the 
cumulative probability to find an isolated region with overdensity 
8c{z) (collapsed region) and mass ^ M (or -y u less than the value 
corresponding to M for the given z). When we apply Eq. l !12t to the 
Gaussian (PS) case, we obtain Pps(7«; z) — 1 — 2F(M, S c \oo, 0). 
Since in such case F is function of j u only, the dependence of 
Ppsi^u', z) on z vanishes (Pps depends on z but through j u ); this 
is the well known fact that for a Gaussian field, the halo multiplicity 
function is universal when expressed in terms of the scaled variable 
v (or 7„). 

By using the Tinker08 halo MF in the first equality of Eq. 
d22t . the cumulative probability function P(7«; z) can be obtained 
numerically. The result is shown in the top left panel of Fig. [8] 
down to the limit mass at each epoch (solid black lines for z = 
0, 1, 2, 3, 4). The corresponding result for the Gaussian (PS) case, 
with S c rescaled by a factor 0.86 in order to fit the high-mass end 
of the Tinker08 MF, is also plotted (red dot-dashed curves). It is 
interesting to note that for the Tinker08 MF, the dependence of 
P(7u) on z is negligible, which suggests that P(~fu) as inferred 
from the simulations deviates only slightly from a universal func- 
tion. In what follows, we will consider that the function P("/u) is 
the same at all redshifts, omitting the explicit argument z. 

In order to make P{"iv) physically meaningful as a cumulative 
probability, we need to know its entire profile, i.e. down to ~y u = 
(M = 0). Then, we proceed first by smoothly extrapolating the 
Tinker08 P("iu) applying a cubic polynomial calculated with the 
conditions to be when 7„ = 0, and fitting the low limit of the 
numerical curve up to the second derivative. The result is shown by 
the dashed black curve in the left upper panel of Fig. [8] The plot 
highlights the difference between the Tinker08 and the PS cases in 
the P{~y u ) shape when 7„ < 0.5. While the latter reaches the ori- 
gin following a straight line, the first one shows a curved shape. The 
physical reason of such difference can be due to the lack of normal- 
ization in the Tinker08 MF; such a lack of normalization can arise 
from the presence of diffuse matter in simulations unaccounted in 
the MF (see below). 

Following, we discuss critical issues that appear when a con- 
frontation between analytical approaches (e.g., ES and CF) and the 
N-body simulations is carried out. Then, taking into account these 
issues, we propose an economical strategy to make compatible for 
the CF, the cumulative probability P(7n) inferred from simulations 
and extrapolated to low masses. In the analytic formalisms, the lin- 




Figure 8. Back integrated cumulative probability P vs. y u (Eq. 1221 for 
the Tinker()8 MF at z = 0, 1, 2, 3, 4 (solid black lines). In each panel the 
virial mass in the MF is rescaled by the factor shown by the label inside. 
The dashed black lines show a cubic extrapolation down to the origin for 
z = 0, fitting the low limit of the cumulative probability up to the second 
derivative. The PS case is displayed in each panel just for visual reference 
(red dot dashed line, see text) 



ear overdense region of mass M is linked to a collapsed halo of the 
same mass (mass conservation; see discussion in §§2.2), and the 
progenitor distributions at each redshift take into account all the 
available mass as part of collapsed halos. The question is whether 
the simulations account for the mass in the same way. 

First, in the simulations the range of halo masses is limited; 
there is no information below the halo mass resolution limit. If one 
extrapolates the fitted halo MF to lower masses and integrate the 
obtained multiplicity function from to infinity, then one finds a 
deficit. This implies that the fitted MF should change its slope (be 
steeper) at smaller masses, and/or that some fraction of the mass is 
actually not in the virialized halos. The latter brings to considera- 
tion a second issue: in the simulations a non-negligible fraction of 
mass is indeed diffuse (not in halos). 

The presence of diffuse matter is due to several reasons: 

(i) A halo is counted only if it contains tens of particles; groups 
with less particles are part of what is considered as diffuse matter. 

(ii) Several authors have shown that a significant fraction of the 
gravitationally bounded particles are actually further away from the 
spherical virial radius (e.g., Prada et al. 2006; Cuesta et al. 2008; 
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Lacerna & Padilla 2011; Anderhalden & Diemand 2011). It is not 
easy to determine the halo radius that contains all the bounded par- 
ticles, and likely this radius presents a large variation from case to 
case, depending on environment, epoch, previous assembly history, 
mass, etc. Let us call ip the average mass fraction of matter grav- 
itationally bounded to halos but not accounted in the conventional 
spherical virial mass. 

(iii) Due to true dynamical processes, e.g., when halos collide, 
some fraction of the particles are ejected from the merged system 
(Wang et al. 2011); let <p be the overall fraction of such ejected 
mass. 

The diffuse mass produced by all of these effects may eventually 
infall on to the growing halos. Therefore, at difference of the ana- 
lytical formalisms, the mass growth of halos in simulations happens 
also in the form of diffuse accretion, which is expected to reduce 
the minor merger rates as compared to the analytical formalisms. 
In the current simulations, at least 30% of the z = halo masses 
came in diffuse accretion (c.f. Genel et al. 201 1; Wang et al. 201 1). 

(iv) It was suggested also a background of diffuse matter, rem- 
nant of the cut-off in the mass power spectrum due to the relativistic 
free-streaming damping (e.g., ~ 10~ 6 M© for the neutralino); at 
very high redshifts, most of this matter due to the cut-off seems to 
be diffuse (Angulo & White 2010), but at the redshifts of signif- 
icant halo mass assembly, the effects mentioned above dominate 
over this. 

The diffuse matter in simulations implies that the mass origi- 
nally linked to a given overdensity is less than that defined inside 
the spherical virial radius. Due to these effects, the conventional 
virial masses result then depleted on average by the fraction if + ip, 

Now, considering the issues discussed above and that in the 
CF the diffuse matter is not taken into account, we need to make as 
compatible as possible the cumulative probability P(7 U ) inferred 
from simulations with the one (correctly normalized from to in- 
finity) to be used in the CF. In the case of item (ii) above, the so- 
lution is simple: a mass rescaling allows us to recover the mass 
inside the collapsed halo proper of the analytical formalism. With 
this correction the analytic approach is justified. The ejected matter 
(item iii) may be considered a simple extension of this approxima- 
tion. After several experiments, we find that an economical way is 
using a strategy based on the following ideas: a) a mass rescaling 
allows to pass from a model (Tinker08) P(~iu) to another P'(7 U ), 
where the complications of the ejected matter and incorrect halo 
mass definition in simulations are avoided, in such a way that the 
masses are now roughly compatible with the ones of the analytic 
case and P'("iu) can be handled with the CF; b) a good empirical 
criterion is to choose P'(7 U ) with a reduced curvature when ex- 
trapolated to low masses in agreement with the shape of Ppsiju)- 
Finally the choice of the mass rescaling factor is justified by the re- 
sult, i.e. the MAHs and MRs being in agreement with simulations. 

Following this strategy, we explore now the effects of a sim- 
ple constant mass rescaling of the Tinker08 halo MF on P'(7 U ). 
The panels of Fig. [8] show the cases of virial mass rescaling by 
factors / = 1.1, 1.2, 1.3, 1.4, and 1.5, respectively. It is note- 
worthy that for / = 1.3 — 1.4, the extrapolation of P'(7 U ) to low 
7u values is close to a straight line and, in general, this probabil- 
ity function approaches to Pps(lu)- Moreover, from the papers 
mentioned above, it can be said that an average halo mass correc- 
tion by tp + ip ~ 30 — 40% (assumed to be mass independent) is 
within the uncertainties. We conclude that this exercise, changing 
the shape of P(7«) obtained from simulations so that it is closer to 
the shape of the analytic (PS) Gaussian probability, implies a mass 



rescaling of / ~ 1.3, in rough agreement with measures of gravi- 
tationally bounded halo mass and diffuse matter considerations in 
simulations. We do not attempt to follow a more rigorous exercise, 
for instance introducing a mass dependence in /, because of the 
large uncertainties involved in the simulation analysis, and because 
at this point we require only an indicative mass correction that en- 
compasses all the complexity implied in passing from the linear 
overdensity to a virialized structure. 



4.2 Mass aggregation histories and merger rates according 
to the mass function from simulations 

The function P("yu) given by Eq. {22} and obtained from the Tin- 
ker08 halo MF is the cumulative unconditional probability to find 
an isolated region with overdensity S c and mass ^ M inside the 
entire universe, being in this case the cosmic -y u — 8 c /\f2oM. 
Now, in order to apply the CF, we need to pass from such uncondi- 
tional probability to a conditional probability. Taking into account 
the generalized PS Ansatz given by Eq. ([5} and according to the 
definition of 7 in Eq. dlOt . the simplest random field compatible 
with the previous reasoning is 



F(M',5'\M,5) 



P(l) 



(23) 



In this way, even if the statistics is not any more Gaussian, Eq. l[5} 
by construction is satisfied and the involved mathematics is simple. 
From Fig. [8] one can see that our inferred P(7) from the Tinker08 
MFs is different from the one corresponding to a Gaussian field 
(see for a similar conclusion Cole et al. 2008; Neistein et al. 2010). 
The probability p(M' , S' \ M, S) dM' to find an isolated region 
with overdensity S' and mass between M' and M' + dM' inside an 
isolated region with overdensity 5 and mass M, according to Eqs. 
{[2]! and l[23]>, is: 



p(M',5'\M,5)dM' 



dP(-y) dy 
dj dM 7 



dM'. 



(24) 



We can apply now the Monte Carlo method to find the halo MAH 
and MRs according to the procedure defined in §§2.2. 

For 2 10 4 Monte Carlo realizations, we calculate the average 
MAHs and MRs by using the cumulative probability P(y) inferred 
from the Tinker08 MF (without any mass rescaling, i.e. / = 1). 
The left panels of Fig. [9] show the average MAH (upper panel, 
solid line) and MRs (lower panel, the same five threshold ratios 
£ as in Fig. [3}. The corresponding average MAH from the Millen- 
nium Simulations (Fakhouri et al. 2010) is shown with the dashed 
line. The predicted MRs are higher at all epochs, but more at lower 
redshifts, than those measured in the simulations by Fakhouri & 
Ma (2008) and Fakhouri et al. (2010; thick ticks in the dN/dz axis 
correspond to the z = MRs as reported by these authors). This 
produces also a too fast mass growth of haloes as seen in the upper 
panel. On the ground of these results, we conclude that the cumula- 
tive probability inferred from the original Tinker08 halo MF leads 
to halo MAHs and MRs in disagreement with those measured in 
numerical simulations. 

In the right panels of Fig. [9] we present the same as in the left 
panels but now for P(7) inferred from the Tinker08 halo MF with 
the virial mass rescaled by the factors / = 1.2, 1.3, and 1.4 (see 
Fig. [8). The solid curves are for / = 1.3 and the left/right (up- 
per/lower) curves are for / = 1.2// = 1.4 in the top (bottom) 
panel. The predicted average halo MAH and MRs are now close 
to those measured in the simulations (certainly within the scatter 
and uncertainties) for the case / = 1.3. It is remarkable that the 
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Figure 9. Left panels: Average MAH (top; solid line) and MRs for different mass-ratio thresholds £ (bottom; color code as in Fig. [5) for a Mq = 10 13 Mq 
present-day halo using the cumulative probability function derived from the Tinker08 MF by Eq. l |23| . The dashed line shows the corresponding mean MAH 
obtained from the Millennium Simulations (Fakhouri et al. 2010). The blue/green thick ticks in the dN/dz-axis show the 2 = MRs for the different £ 
thresholds, as reported in Fakhouri et al. (2010; see Fi^3]. Both the average MAH and MRs do not agree with those measured in the Millennium Simulations 
(see text). Right panels: Same as in left panels but using the cumulative probability function derived from the Tinker08 MF rescaled in mass by / = 
1.2, 1.3, 1.4. The solid lines are for / = 1.3, while the dotted lines for / = 1.2 and / = 1.4 (left and right curves in the top panel, respectively, and upper 
and lower curves for each £ threshold in the bottom panel, respectively). The MAH and MRs agree with those measured in the Millennium Simulation for 
/ ~ 1.3 (see text). 



value / = 1.3 is close to that one which provides a P'("/ u ) with a 
shape similar to the shape of Pps(lu) according to the discussion 
in §§4.1. These two findings are completely independent between 
them. Besides, the value / ~ 1.3 seems to be consistent with mea- 
sures of the diffuse matter component in the simulations (see for 
references §§4.1). Note that the MAH and MRs measured from the 
simulations could change if the virial mass is rescaled; however, we 
do not expect significant changes because in the case of the MAH 
and MR the physics has to do with mass ratios. 

For completeness, we have applied the same procedure to the 
fitting MF given by Sheth & Tormen (1999). We have obtained 
practically the same results as for the Tinker08 MF, which for econ- 
omy we do not repeat here. 

Summarizing, through a mass rescaling that makes compati- 
ble the masses between simulations and the analytic framework, we 
have derived from the halo Tinker08 MF a conditional field proba- 
bility function given by Eq. ( |23t particularly handy in the sense that 
it is compatible with the PS Ansatz (Eq. O, in spite of it deviates 
from the one of a Gaussian density fluctuation field. Starting from 
such probability, our CF allows to calculate the halo MAHs and 
MRs. The obtained average MAHs and MRs are consistent with 



numerical simulations when the virial mass in the input Tinker08 
halo MF is rescaled. 

For economy, we plotted results regarding the average MAH 
and MRs only for one descendant mass, 10 13 Mq. The conclu- 
sions are similar for other masses. Figure [lOjshows the value of / 
required to agree with the Fakhouri et al. (2010) average MAHs at 
the redshift one-half (solid line) and one-tenth (dashed line) of the 
mass for z = descendants of 10 11 , 10 12 , 10 13 , 10 14 , and 10 15 
Mq. As is seen, the mass rescaling factor is roughly the same for 
masses below ~ 10 14 Mq, something in between 1.30 and 1.34, 
in order to agree with the MAHs down to one-tenth of the z = 
descendant mass. For larger masses, this factor should be slightly 
larger. We conclude that, taking into account a mass rescaling of 
1.3, the CF allows for a reasonable good description of the halo 
MFs, and the average MAHs and MRs of halos < 10 14 Mq ob- 
tained in numerical simulations. Probably, a better tuning of this 
description can be attained by varying / with mass, however, this 
is a task dificulted by subtle aspects, for example, the way the aver- 
aging is carried out when calculating the average MAHs and MRs. 
We deserve such an exploration for future works. 

The level of normalization on the Tinker08 MF influences our 
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Figure 10. Mass rescaling factor, /, required to agree with the average 
MAH from Fakhouri et al. (2010) at the redshift one-half (solid line) 
and one-tenth (dashed line) obtained for a wide range of present-day halo 
masses. 



result. Suppose that f rn is a correction factor on the MR The model 
establishes that / = 1.3/ f rn , that is a renormalization on the MF 
produces the same effect of a mass rescaling. If f rn — 1.3, no mass 
rescaling should be necessary at all. A defect in the MF normaliza- 
tion and a defect in the halo mass are related aspects of a same 
problem. If the diffuse matter plays the role that we argued in our 
reasonings above, then the Tinker08 MF does not include all the 
dark matter in the cosmological boxes, and consequently it is not 
normalizable. Tinker08 indeed acknowledge that this is the case. 



5 CONCLUSIONS AND DISCUSSION 

We develop an approach able to connect the overdensities of the 
density fluctuation field with the abundance and assembly history 
of the virialized haloes. Starting from a very general basis, our con- 
ditional formalism (CF) concerns the inventory of the isolated over- 
dense regions given the field conditional probability function of the 
density fluctuation field. This formalism allows us to calculate at 
the same time and in a simple way: the halo mass function MF 
at any z, the halo mass aggregation histories MAHs, and the halo 
merger rate histories MRs. After that, our main goal is to identify 
a strategy that allows us to use the CF for describing optimally the 
numerical simulation results regarding these three halo statistical 
and evolutionary features. In the way of attaining this goal, we ar- 
rived to some important conclusions: 

• For the Gaussian density field (the generalized PS Ansatz, 
Eq. {5J, applies), the predicted average halo MAH and MRs agree 
roughly with those measured in current cosmological (ACDM) 
simulations, but the halo MF, as previously found, is overabundant 
at intermediate masses and deficient at high masses as compared to 
simulations, e.g., with the Tinker08 MF. 

• The introduction of 5 C depending on mass (moving barrier; 
Sheth et al. 2001) instead of the constant one, improves the halo 
MF as compared to simulations, but dramatically spoils the average 
MAH and MRs. Nonetheless, the major problem with the mass- 
dependent 5 C is the excessive sensitivity of the MAH and MRs to 
the redshift step used in the Monte Carlo merger tree construction. 

• A cumulative field unconditional probability function can be 



inferred from the halo MF measured in simulations (Tinker08), be- 
ing this almost the same (universal) at different redshifts when ex- 
pressed in terms of 8 c /o\ its definition has to be complemented by 
adequately extrapolating it to low masses. This inference is par- 
ticularly suitable regarding the shape when the halo virial mass is 
rescaled by / ~ 1.3. Remarkably, this rescaling makes the halo 
mass compatible with the mass used in the analytical formalisms 
having in mind that in simulations not all the mass is in halos. 
On the ground of the CF, the conditional probability function is 
obtained from the unconditional one by making the former com- 
patible with the PS Ansatz, in spite of that the shape of the latter 
deviates from the one corresponding to a Gaussian statistics. Such 
a situation allows us to easily use the CF and our Monte Carlo algo- 
rithm for calculating the halo MAHs and MRs. The obtained aver- 
age MAHs and MRs depend critically on the mass rescaling factor. 
If this factor, assumed constant, is / ~ 1.3, then the agreement with 
the average MAH and MRs in simulations becomes remarkable, at 
least for descendant halos up to ~ 10 14 Mq. It is encouraging that 
the rescaling of the virial mass in the Tinker08 halo MF by a factor 
/ ss 1.3 allows both for a well behaved overall cumulative field 
unconditional probability function and for halo MAHs and MRs in 
agreement with simulations. It is important to remark that the value 
of / depends on the resolution limit of numerical simulations. 

• Dark matter in simulations exhibits a rather complex structure. 
Not all the mass is counted in the halos defined up to the conven- 
tional virial radius, but a diffuse component is present because of 
at least three reasons: a) the mass resolution limit on halo identifi- 
cation, b) the not counting of gravitationally bounded particles that 
are further away the virial radius, and c) the mass ejection generated 
by mergers. The contribution a) has to do with the low mass exten- 
sion of the MF. Due to b) and c) the mass accounted in the virialized 
halo is smaller than in the original overdense region, producing this 
a normalization defect in the MF. The disagreement between the 
simulation and analytic MFs seems to lie on these last two effects 
as well as on a departure from the Gaussian distribution of the den- 
sity fluctuations. The conditional probability function inferred from 
simulations with our approach and its use in the CF allow for an es- 
timate of about 30% for the diffuse mass fraction related to effects 
b) and c), and show that the conditional probability corresponds to 
a density fluctuation field (slightly) deviated from Gaussianity. In 
agreement with the previous argument, the Tinker08 MF does not 
include all the dark matter present in the cosmological boxes and 
consequently it is not normalizable. 

A consistent description of the cosmological simulations re- 
sults regarding the halo MF, MAHs, and MRs through an analytical 
(statistical) formalism has been achieved. The interesting finding is 
that such a description could not be attained for a Gaussian density 
fluctuation field, as well as without taking into account the complex 
distribution of dark matter in simulations. In this sense, our results 
urge for a revision of two questions: 

(i) The correct accounting in simulations of the mass in virial- 
ized halos and of diffuse matter. In the analytical formalisms there 
is an exact equality between the mass of a linear overdensity and 
that of the collapsed halo used for the counting in the MF, while 
in the numerical simulations this is not the case because of effects 
b) and c) above mentioned. Therefore, when confronting the halo 
MF calculated in the analytical formalisms with that one measured 
in simulations, these issues should be taken into account. Our ap- 
proach is a first approximation to this complex problem. Further 
exploration is necessary. 
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(ii) The correct introduction of the statistics in analytical for- 
malisms. The hierarchical structure formation implies that in the 
same density field coexist collapsed structures at small scales with 
regions in the linear regime at large scales. Then, even if the pri- 
mordial statistics is Gaussian, at the time the analytical approach 
is applied the density field at the scales of interest could already 
have deviated from Gaussian initial conditions (e.g., Coles & Jones 
1991), as numerical simulations suggest. 

In any case, the CF and the heuristic approach presented here, 
allow us to overcome some of the difficulties related to these ques- 
tions, and attain a consistent description of the halo MF, MAHs and 
MRs measured in cosmological numerical simulations. 



ACKNOWLEDGMENTS 

We thank the anonymous Referee for his/her comments and sug- 
gestions which helped to improve this paper. C. F. is grateful to 
Dr. Flavio Firmani from the University of Victoria (BC,Ca) for the 
encouraging support to carry out this work. V. A. acknowledges 
CONACyT grant 167332-F for partial funding. 



REFERENCES 

Angulo, R. E„ & White, S. D. M. 2010, MNRAS, 401, 1796 
Anderhalden D., & Diemand J., 2011, MNRAS, 414, 3166 
Avila-Reese V., Firmani C, & Hernandez X., 1998, ApJ, 555, 37 
Avila-Reese V., Firmani C, Klypin A. & Kravtsov A.V., 1999, 

MNRAS, 310, 527 
Avila-Reese V., & Firmani C, 2000, RMxAA, 36, 23 
Bond J.R., Cole S., Efstathiou G., & Kaiser N., 1991, MNRAS, 

379, 440 
Bower R.G., 1991, MNRAS, 248, 332 
Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80 
Carlberg, R. G., & Couchman, H. M. P. 1989, ApJ, 340, 47 
Cole, S., Helly, J., Frenk, C. S., & Parkinson, H. 2008, MNRAS, 

383, 546 
Coles, P., & Jones, B. 1991, MNRAS, 248, 1 
Cuesta A.J., Prada E, Klypin A., Moles M., 2008, MNRAS, 389, 

385 
Fakhouri, O., & Ma, C.-P. 2008, MNRAS, 386, 577 
Fakhouri O., Ma C.-P, & Boylan-Kolchin M., 2010, MNRAS, 

406, 2267 
Firmani C, & Avila-Reese V., 2000, MNRAS, 315, 457 
Firmani C, & Avila-Reese V., 2009, MNRAS, 396, 1675 
Genel, S., Bouche, N., Naab, T, Sternberg, A., & Genzel R., 201 1, 

ApJ, 719, 229 
Gunn J.E., & Gott J.R.I., 1972, ApJ, 176, 1 
Hopkins, P. F. 2012, MNRAS, 423, 2016 
Inoue K.T., & Nagashima M., 2002, ApJ, 174, 9 
Jedamzik K, 1995, ApJ, 505, 37 
Klypin, A., & Rhee, G. 1994, ApJ, 428, 399 
Lacerna I., & PadillaN., 2011, MNRAS, 412, 1283L 
Lacey C, & Cole S., 1993, MNRAS, 262, 627 
Lam, T Y, & Sheth, R. K. 2009, MNRAS, 398, 2143 
Maggiore, M., & Riotto, A. 2010, ApJ, 717, 526 
Mo H.J., van den Bosch E, & White S., 2010, Galaxy Formation 

and Evolution, Cambridge University Press 
Moreno J., Giocoli C, & Sheth R.K., 2008, MNRAS, 391, 1729 
Musso, M., & Paranjape, A. 2012, MNRAS, 420, 369 
Musso, M., & Sheth, R. K. 2012, MNRAS, 423, L102 



Nagashima M., 2001, ApJ, 562, 7 

Navarro J.F., Frenk C.S. & White S.D.M., 1997, ApJ, 490, 493 

Neistein, E., van den Bosch, F. C, & Dekel, A. 2006, MNRAS, 

372, 933 
Neistein, E., Maccio, A. V., & Dekel, A. 2010, MNRAS, 403, 984 
Paranjape, A., Lam, T Y, & Sheth, R. K. 2012, MNRAS, 420, 

1429 
Peacock J.A., & Heavens A.F., 1990, MNRAS, 243, 133 
Prada F. et al., 2006, ApJ, 645, 1001 
Press W.H., & SchechterP, 1974, ApJ, 187, 425 [PS] 
Sheth, R. K, & Lemson, G. 1999, MNRAS, 305, 946 
Sheth R.K., & Tormen G., 1999, MNRAS, 308, 119 
Sheth R.K., Mo H.J. & Tormen G., 2001, MNRAS, 323, 1 
Sheth R.K., & Tormen G., 2002, MNRAS, 329, 61 
Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709 
van den Bosch, F. C. 2002, MNRAS, 331, 98 
Wang, J., Navarro, J. E, Frenk, C. S„ et al. 2011, MNRAS, 413, 

1373 
Wechsler, R. H., Bullock, J. S., Primack, J. R., Kravtsov, A. V., & 

Dekel, A. 2002, ApJ, 568, 52 
Yano T, Nagashima M., & Gouda N., 1996, ApJ, 466, 1 
Zentner A.R., 2006, |astro-ph/061 1454 1 
Zhang J., Fakhouri O., & Ma C.P., 2008, MNRAS, 389, 1521 



